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1.0  SUMMARY 


This  report  addresses  a  series  of  simulations  of  turbulent,  multi-species  flow  fields 
surrounding  a  sphere.  The  purpose  of  these  simulations  is  to  validate  our  capability  for  using  the 
Large  Eddy  Simulation  with  Linear  Eddy  Modeling  in  3  Dimensions  (LESLIE3D)  multiphase 
physics  computer  program.  A  high  Reynolds  number  flow  field  is  computed  at  Mach  2.0,  and  a 
low  Reynolds  number  flow  field  is  calculated  at  Mach  0.1.  The  former  flow  field  is  strongly 
shocked  while  the  latter  flow  field  is  shock  wave  free.  In  each  test  calculation,  the  physics  of  the 
flow  field  is  examined,  and  the  drag  coefficients  are  calculated  via  post-processing  for  comparison 
with  data.  At  the  same  time,  drag  post-processor  algorithms  are  validated. 
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2.0  INTRODUCTION 


During  the  past  three  decades,  research  in  the  field  of  computational  fluid  dynamics  (CFD) 
has  shown  remarkable  progress.  CFD  computer  codes  are  now  widely  applied  in  the  commercial 
world  for  aircraft  design  with  little  requirement  for  wind  tunnel  testing.  A  wide  range  of  flight 
conditions  for  subsonic,  transonic,  supersonic  and  hypersonic  flight  can  now  be  addressed  by 
continuum  based  CFD  codes  at  nominal  cost.  The  inclusion  of  more  complex  equations  of  state 
and  finite  rate  chemical  reaction  algorithms  allow  the  physics  of  propulsion  systems  to  be  solved 
such  as  the  burning  of  fuel  in  gas  turbine  combustors.  Intricate  multiphase  physics  equations 
couple  the  behavior  of  gas  phase  CFD  algorithms  with  dispersed  fields  of  solid  particles  and/or 
liquid  droplets.  This  innovation,  associated  with  more  recent  efforts,  allows  the  numerical  analysis 
of  the  environment  within  burning  solid  rocket  motors.  Because  of  its  influence  on  chemistry  and 
on  dispersed  particle/droplet  fields,  turbulence  becomes  a  dominant  factor  in  the  solution  of 
engineering  problems.  Within  the  past  decade,  Large  Eddy  Simulation  (LES)  has  evolved  to 
become  a  major  force  in  solving  many  types  of  multiphase  problems.  Modem  LES  advances 
beyond  the  older  Reynolds  Averaged  Navier-Stokes  (RANS)  turbulence  models  by  exploiting 
mathematical  similarity  between  the  resolved  Leonard  stress  and  the  modeled  subgrid  stress. 
Through  the  use  of  proper  closure  terms,  modern  LES  algorithms  can  now  capture  the  correct  rates 
for  finite  rate  chemical  reactions.  All  of  these  capabilities  exist  within  LESLIE3D,  the  Large  Eddy 
Simulation  with  Linear  Eddy  model  in  3  Dimensions,  a  computer  code  developed  by  Professor 
Suresh  Menon  at  the  Georgia  Institute  of  Technology. [1] 

Guided  by  the  technological  needs  of  the  Munitions  Directorate,  over  the  past  fourteen 
years  LESLIE3D  has  evolved  from  being  designed  for  subsonic  combustion  applications  to  a 
versatile  code  suited  for  all  flight  regimes  with  dense  multiphase  fields.  Its  crowning  achievement 
is  found  in  its  shock  turbulence  algorithms.  In  earlier  days,  shock  capturing  algorithms  held  a 
distinct  conflict  of  interest  with  turbulence  simulation  algorithms.  Shock  capturing  methods  are 
overly  dissipative  tending  to  wash  subtle  turbulent  fluctuations  out  of  the  numerical  flow  field.  On 
the  other  hand,  turbulence  simulation  schemes  often  apply  centered  difference  stencils.  Centered 
stencils  do  not  respect  the  directionality  of  acoustic  waves  in  the  flow  field  causing  destructive 
oscillations  near  shock  waves.  Careful  research  culminating  in  2005  produced  the  hybrid  scheme 
employed  by  LESLIE3D  to  combine  a  set  of  low  dissipation  shock  capturing  schemes  with  the 
highly  accurate  Locally  Dynamic  subgrid  Kinetic  energy  Model  (LDKM)  developed  by  Professor 
Suresh  Menon.  These  algorithms  are  core  capabilities  installed  within  the  parallel,  multi-block 
version  of  LESLIE3D  applied  for  the  problems  discussed  here. 

An  issue  that  routinely  confronts  computational  physics  computer  codes  is  that  of  code 
validation.  Without  getting  into  details,  validation  is  intended  to  ensure  that  the  computer  solution 
generated  for  a  particular  problem  renders  results  that  compare  favorably  either  with  experimental 
results  or  with  archived  or  exact  solutions.  For  the  standpoint  of  the  computational  physicist, 
physics  computer  codes  should  be  developed  based  on  “unit  experiments”  that  elucidate  the  core 
physics  behind  the  problem  in  question.  These  experiments  are  highly  controlled  and  above  all, 
have  repeatable  results  that  form  the  foundation  of  a  validation  dataset.  The  impression  of 
validation  held  by  the  project  engineer  (or  the  customer)  is  often  different.  The  customer  often 
views  validation  as  involving  the  comparison  of  code  results  for  a  much  larger,  less  well  controlled 
test  configuration.  In  many  cases,  the  type  of  test  performed  is  not  truly  repeatable.  The  statistical 
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variation  for  measured  parameters  may  be  much  larger.  As  a  result,  comparisons  between  test  and 
computer  code  results  may  differ  to  a  greater  degree.  In  many  cases  the  differences  result  from 
unknown  or  imprecisely  known  aspects  of  the  test.  For  example,  the  precise  location  of  a  solid 
wall  or  say,  an  accurate  representation  of  the  surface  roughness  of  a  wall  may  not  be  known. 
Secondly,  the  same  parameter  may  change  between  different  test  realizations.  Addressing  this  type 
of  uncertainty  goes  beyond  the  basic  computational  physics  equations.  The  field  of  uncertainty 
quantification  (UQ)  is  required  to  address  these  issues.  The  interested  reader  is  referred  to  [2]  for 
a  discussion  of  this  field. 

At  first  glance,  validation  certainly  encompasses  the  performance  of  the  computational 
physics  computer  code,  but  it  also  must  address  the  post-processing  algorithms  used  to  transform 
code  results  into  usable  engineering  parameters.  As  an  example,  many  CFD  codes  produce  files 
containing  say,  pressure,  temperature  and  velocity  components  at  each  point  in  the  flow  field. 
These  are  critical  physics  properties  of  the  flow  field,  but  they  are,  in  many  cases,  not  directly 
useful  to  the  engineer.  Instead,  the  engineer  requires  knowledge  of  the  force  components  (lift  and 
drag)  or  aerodynamic  moments.  Also,  it  may  be  desirable  to  know  the  rotational  character  of  the 
flow.  In  each  case,  the  data  provided  by  the  computational  physics  computer  program  must  be 
post-processed  to  extract  the  required  data.  This  process  involves  algorithms  that  also  require 
validation. 

Our  outstanding  project  needs  involve  the  simulation  of  high  speed,  turbulent  flow  around 
solid  bodies  where  flow  separation  may  occur.  It  is  important  to  compute  forces  on  this  type  of 
configuration,  so  we  have  selected  a  sphere  as  the  solid  body,  and  our  task  is  to  compute  the  time 
varying  forces  exerted  on  it  by  the  flow  field.  By  comparing  the  average  drag  coefficient  exhibited 
by  the  sphere  against  archival  data,  we  obtain  a  validation  point  for  the  post-processing  algorithms 
and  for  the  user's  LESLIE3D  input.  Documenting  the  results  of  this  validation  effort  is  the  goal  of 
this  report. 
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3.0  METHODS,  ASSUMPTIONS  AND  PROCEDURES 


3.1  Navier-Stokes  Equations 

An  assessment  of  the  physics  of  high  speed  flow  fields  under  normal  atmospheric 
conditions  (excluding  free  molecular  flow  in  rarefied  air)  for  common  finite  length  scales  implies 
the  existence  of  high  Reynolds  number  flow.  In  a  simple  interpretation,  the  Reynolds  number,  a 
dimensionless  quantity,  represents  the  ratio  of  inertial  forces  to  viscous  forces  in  the  flow  field.  [3] 
The  Reynolds  number  is  defined  as 


Re  =  ^ 


(1) 


where  p  and  V  are  the  local  flow  density  and  speed.  L  is  the  characteristic  length  scale,  and  p  is 
the  dynamic  viscosity.  Low  Reynolds  numbers  are  indicative  of  flows  where  viscous  effects  are 
dominant;  turbulence  is  unlikely.  In  high  Reynolds  number  flows,  inertial  effects  are  important 
and  overwhelm  viscous  effects.  These  flow  fields  tend  to  incur  instability  that  leads  to  the 
formation  of  turbulence.  All  real  flow  fields  are  governed  by  the  conservation  of  mass,  momentum, 
energy  and  chemical  elements.  In  the  continuum  limit,  these  properties  are  conserved  by  the 
Navier-Stokes  equations  [4],  i.e., 


f  +  v(Pv)  =  0 


(2) 


is  the  conservation  of  mass  (or  continuity  equation)  where  p  is  the  density  of  the  gas  mixture  and 
V  is  the  flow  velocity  vector.  The  time  coordinate  is  1.  The  momentum  equation  may  be  written  as 

+  V-(pVV)  +  VP  =  V-t  (3) 

In  equation  (3),  the  construct  VV  is  the  tensor  (or  dyad)  product  of  the  velocity  vector  with  itself 
resulting  in  the  second  order  advection  tensor.  The  thermodynamic  pressure  is  denoted  by  P.  The 
final  term  in  this  equation  is  t,  the  deviatoric  stress  tensor.  It  is  in  this  term  that  the  viscous  effects 
are  included.  This  tensor  may  be  written  in  Cartesian  components  as 


r 


-i,j 


auA 

dxj 


(4) 


where  the  velocity  vector  is  denoted  as  V  —  (it1;  u2,  u3)T .  The  parameter  A  is  denoted  as  the  bulk 
viscosity  coefficient.  Using  Stake's  hypothesis,  this  coefficient  is  usually  expressed  as 


A  = 


(5) 


Equation  (3)  is  an  expression  of  Newton's  second  law  for  the  flow  field  in  the  absence  of  body 
forces.  Note  that  equation  (3)  is  vector  valued.  It  actually  represents  three  distinct  scalar  equations. 
The  conservation  of  total  energy  is  enforced  by  the  energy  equation,  i.e., 
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d<?i 

dxi 


(6) 


d(pE)  d(pujE)  djujP )  _  d{xjj) 

dt  dxj  dxj  dxi 


where  E  is  the  total  energy  per  unit  mass  given  by 

E  —  e  + 


(7) 


Note  that  the  summation  convention  is  in  effect  over  repeated  indices.  Also,  the  internal  energy 
per  unit  mass  is  denoted  as 


=  e°  +  Sr0CVin(f)df  (8) 

For  a  single  calorically  perfect  gas,  it  is  advantageous  to  use  the  form 

e  =  CVT  (9) 

For  species  n,  the  constant  volume  specific  heat  is  denoted  Cv  n  while  T is  the  absolute  temperature. 
The  final  term  in  equation  (6)  represents  the  heat  flux  vector  defined  within  the  fluid;  specifically, 

qi  =  +  p  En=i  Yn  KVi,n  (10) 

where  k  is  the  thermal  conductivity  of  the  fluid,  and  hn is  the  sensible  enthalpy  of  species  n 
computed  as 


K  =  K  +  fioCPjn(f)df  (ii) 

The  system  comprised  of  equations  (2),  (3)  and  (6)  is  the  common  system  of  conservative  Navier- 
Stokes  equation,  but  to  achieve  mathematical  closure,  the  system  must  include  an  equation  of  state 
such  as  that  for  a  perfect  gas,  i.e., 

P  =  pRT  (12) 

R  is  the  perfect  gas  constant;  that  is 


R 


Eunlv 

MW 


(13) 


In  this  case,  R  is  written  for  a  single  species  with  molecular  weight  MW.  For  a  mixture  of  N  gaseous 
species,  the  perfect  gas  constant  is  written  as 


Rr 


=  R, 


yiV 
En=l 


MWr, 


(14) 


where  Yn  is  the  mass  fraction  of  the  nth  species,  and  Runiv  is  the  universal  gas  constant.  The 


problems  considered  in  this  report  involve  more  than  one  gaseous  species.  LESLIE3D  has  been 
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generalized  to  include  multiple  species,  so  the  system  of  governing  equations  is  augmented  by  a 
species  equation  that  is  written  as 


(15) 


where  Vi  n  is  the  diffusion  velocity  defined  by 


(16) 


Dn is  the  diffusion  coefficient  for  species  n.  [5] 


In  equation  (4),  viscosity,  a  key  property  for  viscous  flow  fields,  is  computed  as 


(17) 


a  relationship  known  as  Sutherland's  law  where  /i0  and  T0  are  the  reference  viscosity  and 
temperature,  respectively,  and  S  is  a  constant.  These  parameters  depend  on  the  gaseous  substance 
in  question.  To  attain  highly  accurate  results  for  mixtures  of  gases,  a  viscous  transport  model  may 
be  needed. 

3.2  Filtered  Navier-Stokes  Equations  for  Large  Eddy  Simulation 

Fundamentally,  the  Navier-Stokes  equations  are  based  upon  the  continuum  hypothesis 
where  it  is  assumed  that  mass  is  modeled  as  a  mathematical  continuum,  not  as  a  vast  collection  of 
tiny  particles. [6]  This  assumption  works  quite  well  as  long  as  the  mean  free  path  for  molecular 
motion  remains  far  smaller  than  the  characteristic  dimension  of  the  immersed  body.  This 
assumption  holds  extremely  well  for  the  application  described  in  this  report.  One  may  ask  that, 
given  the  viability  of  these  equations,  why  not  solve  them  directly  for  the  fluid  properties  of 
interest?  Why  resort  to  the  added  complexity  of  large  eddy  simulation?  The  answer  to  these 
questions  is  reliant  on  required  grid  resolution  and  its  impact  on  computer  memory.  The  reader 
may  recall  that  the  direct  numerical  solution  of  the  Navier-Stokes  equation  is  referred  to  as  Direct 
Numerical  Simulation  (DNS).  For  DNS,  if  one  knows  the  initial  conditions  well  enough,  the 
Navier-Stokes  equations  could,  in  theory,  be  solved  to  produce  a  highly  accurate  solution.  The 
difficulty  associated  with  this  idea  lies  in  grid  resolution.  Turbulence  has  a  range  of  scales  that 
must  be  resolved  to  grant  an  accurate  solution. [7]  Accordingly  ,the  grid  must  be  resolved  well 
enough  to  capture  the  smallest  scales.  This  grid  resolution  scale  as  the  Reynolds  number  raised  to 
the  9/4  power.  [8]  That  is  to  say,  a  simulation  where  Re  =  8  million  requires  a  grid  possessing  about 
3.4  x  1015  points.  This  figure  still  remains  well  beyond  the  capability  of  current  computer  designs. 
So  DNS  is  confined  to  very  small  Reynolds  numbers.  Since  DNS  is  not  a  viable  option,  we  must 
exploit  the  separation  between  the  scales  of  fluid  and  filter  the  Navier-Stokes  equations  in  three- 
space.  In  doing  so,  we  can  separate  the  large  scales  of  motion  (for  resolution  via  the  governing 
equations)  from  the  small  scales  of  motion.  Fluid  activity  at  the  smaller  scales  is  modeled;  thus, 
the  large  eddy  simulation  requires  substantially  less  computational  resources  than  does  DNS.  The 
mathematical  development  of  the  filtered  equations  is  lengthy,  so  only  the  results  are  stated  here. 
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For  an  in-depth  discussion  of  the  derivations,  the  interested  reader  is  referred  to  [5]. 

The  filtered  mass  conservation  equation  is  written  as 

iR  +  «2S>  =  0  (18) 

dt  dxi  v 

where  the  overbar  indicates  the  spatial  (top  hat,  in  this  case)  filter.  The  tilde  notation  indicates  the 
use  of  mass  (Favre)  averaging,  i.e., 


Ui 


pui 

p 


(19) 


The  filtered  momentum  equation  is  derived  by  a  similar  process;  this  equation  contains  the 
modeled  term  rsf s  and  is  written  as 


d(p  u{) 
dt 


+  4r.  [pUiUj  +  PSij  +  T i9jS  ~  T ij]  =  0 


More  specifically,  the  subgrid  stress  tensor  is  given  by 


(20) 


f 


\ 


Tu  =P 


//,  U  ,  —  Uj  li  j 


(21) 


V  J 


Due  to  the  presence  of  a  correlation  term,  the  subgrid  stress  (21)  must  be  modeled.  The  filtered 
total  energy  equation  may  be  written  as 


d(pE) 

dt 


Ui^ji 


+  H-as  +  u-9S ]  =  0 


(22) 


where 


q,  =  ~k{T)  |  +  pS.  i  Ynhn(T )  VLn  +  qfs  (23) 

The  modeled  terms  in  this  case  are  the  subgrid  total  enthalpy  flux,  the  subgrid  viscous  work  and 
the  subgrid  heat  flux  qJ5S.[l]  Respectively,  these  terms  are  mathematically  defined  as 


nr=p 


Eu i  -Eui 


+ 


PUj  -  Pu , 


(24) 


7 

Distribution  A 


The  filtered  total  energy  per  unit  mass  is  given  by 


E  =  e  +  \ukuk  +  ksgs 

where  the  subgrid  kinetic  energy  is  defined  by 


(26) 


k 


SgS  _  J_ 
2 


M'k^'k 


(27) 


We  can  write  the  filtered  species  equation  as  shown  below. 


d(pYn)  | 

dt 


U:  + 


n  =  l,...,N 


(28) 


The  two  subgrid  closure  terms  existing  in  equation  (26)  are  mathematically  defined  as 


y  sgs 
1  i,n 


-P 


U;  Y„ 


(29) 


the  subgrid  convective  species  flux  and 


0%  =P 


V.  Y  -V.,Y 

i,n  n  i,k  n 


(30) 


the  subgrid  diffusive  species  flux.[l]  As  is  the  gas  with  the  other  governing  equations,  the  multi¬ 
species  perfect  gas  equation  of  state  also  requires  filtering.  The  result  is 

P  =  p(Rf  +  RmJ™)  (31) 

This  equation  of  state  contains  Tsgs ,  the  temperature-species  correlation  given  mathematically  as 
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YT-YT 


(32) 


The  internal  energy  must  also  be  filtered.  It  can  be  shown  for  calorically  perfect  gases  that 

N 

‘  +yM„)  (33) 

n—\ 

where  Ah' fn  is  the  standard  heat  of  formation  at  T°,  a  reference  temperature.  Specifically, 

A h'f,n  =  Ahf  n  -  CPnT°  (34) 


3.3  LES  Modeling  Considerations  for  the  Filtered  Governing  Equations 

A  significant  investment  in  research  time  has  been  made  to  adapt  the  LES  framework 
briefly  discussed  above  for  fully  compressible  flow  fields.  As  it  happens,  a  key  aspect  involved  in 
formulating  a  compressible  model  entails  altering  the  evolution  equation  for  the  subgrid  kinetic 
energy.  Although  this  equation  acts,  in  some  sense,  like  a  governing  equation,  it  is,  in  effect,  a 
model,  so  it  presented  in  this  section  of  the  report.  An  earlier  form  [9]  of  the  subgrid  kinetic  energy 
equation  is  written  as 


dpksgs 

dt 


+  -?-(pUikS3S)  =  PS3S  -  Ss3s  +  /-(p 

C /  A  j  U  A  j  \ 


vt  dksas 
Prt  dxt  . 


where  the  production  and  dissipation  terms  are,  as  one  may  suspect,  modeled  as 


psgs 


sgs  du± 
iJ  dxj 


£  — 


(35) 


(36) 

(37) 


The  parameter  vtis  the  turbulent  eddy  viscosity,  and  Prt  is  the  turbulent  Prandtl  number.  The 
width  of  the  filter  is  given  by  A,  and  C£is  a  parameter  that  is  automatically  adjustable. 


As  one  may  suspect,  each  term  in  this  equation  has  a  specific  meaning.  The  first  and  second  terms 
on  the  left  side  of  the  equation  are  the  local  and  convective  rates,  respectively.  On  the  right  side  of 
the  equation,  we  have  the  production,  dissipation  and  diffusion  terms,  respectively,  for  subgrid 
kinetic  energy.  Equation  (33)  is  suited  moreover  for  lower  speed  flows  where  compressibility 
effects  are  small.  Modifications  are  required  for  the  inclusion  of  major  compressibility.  A 
modified  form  that  does  include  Mach  number  effects  may  be  written  as 
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The  terms  on  the  right  side  of  (36)  are  the  diffusion,  pressure-dilatation  correlation,  production 
and  dissipation,  respectively,  for  subgrid  kinetic  energy.  The  diffusion,  pres  sure -dilatation 
correlation  and  dissipation  cannot  be  directly  computed  from  resolved  scale  quantities,  so  these 
terms  require  modeling  based,  in  part,  upon  their  mathematical  definitions.  In  the  interest  of 
brevity,  only  the  final  subgrid  kinetic  energy  equation  is  presented  here. [5] 


d  d 

—  pksgs  _j_  — (puj/cs5s)  = 

a  +  OXj 


dt 


dx; 


( pvt  +  M) 


dksgs  pvtR  df 


dXi 


+ 


1  +  «*(«?")’(*£=)' 


'  sgsduj^ 

Si  dxt  +  PLe 


Prt  dxi 

(ksgs^3/2 


(39) 


where  vt  is  the  turbulent  eddy  viscosity;  p  is  the  dynamic  viscosity,  and  R  is  the  mass-averaged 
mixture  perfect  gas  constant.  Mstas  is  the  turbulent  Mach  number  based  upon  ksgs .  Also,  we  have 
that 


where 


Stj 


1  ( dui 

2  \dxj 


+ 


aiq\ 

dxi) 


(40) 


(41) 


A  major  difference  existing  in  the  new  subgrid  kinetic  energy  equation  (38)  is  that  the  pressure- 
dilatation  correlation,  as  shown  as  the  second  line  in  the  equation,  is  expressed  as  a  Mach  number 
dependent  quantity;  the  parameter  cr^is  a  model  coefficient.  A  new  expression  for  the  dissipation 
of  turbulent  subgrid  kinetic  energy,  Dksgs ,  is  utilized,  i.e., 


Dksgs 


l bc£(ksas )3/2 

A 


(42) 


As  described  above,  the  subgrid  kinetic  energy  is  an  important  property  in  modeling 
turbulence  in  the  flow  field.  Directly  dependent  upon  this  property  is  the  subgrid  stress  tensor.  It 
is  modeled  as  follows. 


<f  =  -2 pvt(Su  -  jSw«y)  +  (43) 

With  subgrid  kinetic  energy  provided  by  the  solution  to  (38),  (42)  becomes  mathematically  closed. 
Basic  model  coefficients  may  be  obtained  for  the  turbulent  eddy  viscosity  and  dissipation.  Note 
that 
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(44) 


vt  =  CvVF^d 


and 


e 


C£ 


{U.sgs^/2 

A 


(45) 


An  examination  of  turbulent  spectra  [10]  shows  that  good  values  for  Cv  and  Ce  can  be  estimated 
as 


Cv  =  0.067  ;  C£  =  0.916  (46) 

Other  values  are  possible  for  these  coefficients,  and  dynamic  evaluation  of  these  coefficients  is  a 
validated  option  in  LESLIE3D.  Since  this  variant  is  not  applied  in  the  present  work,  it  is  not 
discussed  here. 


There  are  other  modeling  terms  that  arise  from  the  filtering  process.  A  brief  discussion  of 
these  terms  is  warranted.  Recall  the  filtered  total  energy  equation  (21).  The  closure  terms  that  arise 
in  this  equation  are  the  subgrid  total  enthalpy  flux  H-gs( 22)  and  the  subgrid  viscous  work  <j-gs(23); 
accordingly,  they  are  modeled  [5]  as  a  sum,  i.e., 


H?gs  +  0;9S 


-  (pvt  +  P) 


dks9s 


dxi 


+  u,xS9S 

Prt  dxi  J  l’J 


(47) 


All  of  the  terms  on  the  right  side  of  (45)  are  either  already  modeled,  known  parameters  or 
properties  extracted  from  the  resolved  field.  The  first  two  terms  on  the  right  side  are  developed  by 
a  process  known  as  gradient  modeling  since  each  involves  a  gradient  computed  from  another  flow 
property.  A  similar  form  exists  for  the  subgrid  convective  species  flux  F^sin  (28).  That  is 


Ysas 

Ii,n 


(48) 


In  the  above  equations,  one  notes  the  presence  of  the  turbulent  Schmidt  number  Sct,  the  ratio  of 
the  turbulent  viscous  diffusion  rate  to  the  molecular  diffusion  rate.  This  dimensionless  number  is 
taken  as  unity.  [1]  The  subgrid  species  diffusive  flux  term  S^in  (29)  is  regarded  as  small  in  high 
Reynolds  number  flow  fields,  so  it  is  neglected.  The  temperature-species  correlation  term  Tsgs 
(30)  is  also  believed  to  be  very  small  in  magnitude,  so  it  too  is  neglected  as  is  the  subgrid  heat  flux 
qsgs  introduced  in  (22). 

The  equations  shown  in  Section  2  represent  a  translation  of  physical  laws  (conservation 
equations,  equation  of  state)  into  the  language  of  mathematics.  This  process  serves  to  formulate  a 
mathematical  model  for  fluid  flow,  but  it  is  not  without  assumptions.  This  system  of  equations  is 
built  upon  the  continuum  hypothesis,  and  although  it  is  somewhat  tacit,  a  Cartesian  coordinate 
system  resides  under  the  mathematics.  It  is  also  worth  mentioning  that  the  perfect  gas  equation  of 
state  is,  in  fact,  a  model,  an  interpretation  of  reality  not  entirely  unlike  a  drawing  made  by  an 
illustrator.  Regarding  the  scenarios  for  which  they  are  designed,  these  equations  are  sound  and 
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form  a  reliable  and  quantifiable  mathematical  model.  Still,  one  should  never  lose  sight  of  the  fact 
that  this  system  of  equations  is  only  a  sketch  of  reality,  albeit  a  detailed  one.  In  Section  3,  we  learn 
that  it  is  not  currently  possible  for  us  to  extract  the  level  of  detail  inherent  in  this  mathematical 
model.  A  primary  inhibiter  is  the  lack  of  computer  memory.  A  more  surreptitious  difficulty  is  that 
at  the  same  scales,  we  may  be  unable  to  accurately  describe  both  boundary  and  initial  conditions. 
For  this  reason,  we  back  away  from  the  fine  scale  and  take  a  more  macroscale  view  of  the  problem. 
That  is  to  say,  we  apply  filters  to  the  mathematical  flow  field  and  extract  “averaged”  or  “filtered” 
flow  field  properties  in  space.  This  process  lacks  neatness  in  that  it  introduces  closure  terms  that 
must  be  modeled  via  a  set  of  auxiliary  equations  that  sit  aside  from  the  conservation  laws  requiring, 
in  some  cases,  the  use  of  empirical  data  and  curve  fitting.  In  the  mathematical  sense,  closure  terms 
add  variables  to  the  system.  For  that  reason,  we  must  add  equations  to  the  system  to  affect  closure 
and  admit  a  solution.  Section  4  introduces  specific  closure  term  models  including  the  time  and 
space  dependent  model  for  subgrid  kinetic  energy,  a  property  used  to  elucidate  the  subgrid  stress 
tensor.  In  fact,  subgrid  kinetic  energy  is  governed  by  an  evolution  equation.  It  is  interesting  to 
realize  that  although  this  equation  has  a  production  term,  it  is  not  capable  of  predicting  the 
development  of  turbulence  where  none  exists  as  an  initial  value.  This  equation  is  “sourceless”. 
One  must  prescribe  a  turbulent  field,  one  with  non-zero  subgrid  kinetic  energy  and  non-zero 
velocity  components,  a  priori.  A  non-turbulent  flow  field  simply  evolves  with  zero  turbulence. 
Then  one  may  ask,  of  what  good  is  the  subgrid  kinetic  energy  equation  if  it  cannot  predict  the 
creation  of  turbulence?  The  answer  to  this  question  is  provided  through  an  understanding  of  the 
evolution  equation.  The  evolution  equation  describes  how  a  property  evolves  from  a  presumed 
initial  value.  Suppose  that  a  non-physically  weak  field  of  turbulence  is  prescribed  at  the  initial 
condition.  If  so,  the  evolution  equation  alters  this  field  based  upon  other  physical  inputs  and 
strengthens  the  turbulent  field.  On  the  other  hand,  suppose  that  a  turbulent  field  is  prescribed  where 
little  or  no  turbulence  really  exists.  In  this  case,  based  upon  other  physical  inputs,  in  time  the 
evolution  equation  depletes  away  the  turbulent  field  as  though  none  ever  existed;  the  dissipation 
term  dominates  the  equation.  The  subgrid  kinetic  energy  equation  endeavors  to  ensure  that  the 
evolving  turbulent  field  is  physically  realizable. 

3.4  Numerical  Algorithms 

In  the  preceding  three  sections  of  this  report,  we  introduced  and  discussed  the  Navier- 
Stokes  equations  for  multi- species,  non-reacting  flow  fields.  Then  spatial  filtering  is  applied  to  the 
equations,  and  the  resulting  closure  terms  are  modeled.  This  development  is  largely  mathematical; 
we  have  not  discussed  any  of  the  numerical  methods  associated  with  solving  these  equations.  In 
this  section  of  the  report,  we  present  a  brief  discussion  of  the  core  solver  methodology;  the  details 
are  omitted  because  a  full  description  of  LESLIE3D’s  numerical  algorithms  is  very  lengthy  and 
covers  many  years  of  research.  For  these  details,  the  interested  reader  is  referred  to  References  [5] 
and  [10]. 

It  is  worthwhile  to  sketch  the  basic  layout  of  LESLIE3D.  The  version  used  to  support  the 
present  work  is  referred  to  as  the  “SVN  code”,  a  loose  designation  named  for  the  source  code 
control  method  used  for  its  development.  The  SVN  code  is  admittedly  an  older  computer  code  that 
is  used  predominantly  for  research,  and  it  contains  the  primary  computational  physics  algorithms 
for  solving  a  wide  array  of  problems.  A  principal  capability  of  LESLIE3D  is  that  it  is  a  multi-block 
computer  program.  The  computational  domain  is  broken  apart  into  “blocks”  of  grid.  As  a  result, 
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LESLIE3D  can  solve  problems  on  irregular  geometries  such  as  winged  bodies  like  aircraft,  or  on 
the  spherical  shell  used  for  problems  addressed  later  in  this  report.  Through  the  Message  Passing 
Interface  (MPI),  the  data  associated  with  these  blocks  can  communicate  with  more  blocks  being 
assigned  to  a  processor  or  core.  It  follows  that  LESLIE3D  is  also  fully  parallelized.  It  ports  and 
operates  in  parallel  to  a  wide  variety  of  computer  systems  including  LINUX  and  CRAY  clusters. 
Problems  of  local  interest  have  been  run  on  over  4,000  cores. 

A  structured  grid  must  be  generated  for  each  block.  For  three-dimensional  problems,  this 
requirement  implies  that  there  is  a  local  (i,  j,  k )  index  system  assigned  to  the  block,  so  each  grid 
point  has  its  own  ordered  index  triplet  (i,  j,  k).  One  can  “navigate”  through  the  grid  in  simple 
counting  order  where 

i  =1,  2,  ...,  IMAX;  j  =  1,  2,  ...,  JMAX,  andk=  1,  2,  ...,  KM  AX. 

The  index  maxima  IMAX,  JMAX  and  KMAX  can,  with  some  restrictions,  differ  for  each  block. 
Moreover,  it  follows  that  each  block  has  six  sides  forming  a  hexahedron  in  the  computational 
plane.  Specifically,  a  given  block  has  a  side  where  i  =  1  and  a  side  where  i  =IMAX.  The  other 
sides  of  the  block  correspond  to  indices  j  and  k.  It  is  on  each  of  these  block  sides  (or  interfaces) 
that  either  boundary  conditions  are  enforced  and/or  parallel  communication  is  established.  A 
restriction  on  maximum  block  indices  arises  at  the  block  sides.  LESLIE3D  presumes  that  block 
grids  are  constructed  point-on-point.  This  statement  implies  that  if  two  blocks  communicate  on  a 
particular  interface,  then  both  blocks  possess  the  same  points  on  the  interface.  As  a  result,  both 
blocks  must  share  the  same  maximum  lateral  indices  along  the  side  shared  by  the  blocks.  That 
having  been  written,  the  “direction”  of  say,  index  i  in  one  block  need  not  be  the  same  as  an 
adjoining  block.  This  property  also  applies  to  the  other  indices.  Another  restriction  is  that  the  index 
systems  in  each  block  must  be  right-handed  in  the  physical  space. 

One  may  correctly  surmise  that  LESLIE3D  is  a  structured  finite  volume  code.  The 
hexahedral  topology  of  a  computational  block  extends  to  its  interior,  so  the  block  is  comprised  of 
hexahedral  cells.  Since  the  grid  is  structured,  if  a  block  has  maximum  grid  indices  IMAX,  JMAX 
and  KMAX,  then  the  cell  indices  will  have  maxima  IMAX-1,  JMAX-1  and  KMAX-1.  It  follows 
that  the  total  number  of  cells  for  the  block  is  (IMAX- 1  )x(JM AX- 1  )x(KM AX- 1 ) .  Neither  the  blocks 
nor  the  finite  volumes  need  to  be  brick-shaped.  LESLIE3D  easily  allows  arbitrary  geometries;  a 
Cartesian  grid  system  is  not  required.  LESLIE3D  computes  cell  normal  vectors  and  metrics  as 
required  for  both  the  advective  and  viscous  terms.  A  general  curvilinear  coordinate  transformation 
is  inherent  in  the  coding.  A  strength  of  the  finite  volume  method  is  noted  when  examining  the 
governing  equations.  Although  significant  manipulation  is  required  to  achieve  it,  the  system  of 
governing  equations  can  be  written  in  the  vector  form 

^  +  V-F  =  0  (49) 

Ot 

In  this  expression,  Qis  the  vector  of  conserved  variables  i.e., 

Q  =  { p ,  pu,  pv,  pw,  pE,  pksgs,  pYlt ... ,  pYN}T  (50) 
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The  construct  F  is  a  vector  of  flux  vectors  containing  all  of  the  advective  and  viscous  terms 
found  in  the  governing  equations.  It  may  be  expressed  as 


F  = 


iFxJFy,kFg 


(51) 


where  t  f  and  fc  are  the  Cartesian  unit  vectors.  On  the  other  hand,  vectors  such  as  Q,  Fx,  Fy  and  Fz 
are  not  geometric  vectors  in  three-space  (9f3);  rather,  they  are  vectors  in  the  sense  that  they  are 
linear  arrays  of  length  6+N.  To  see  the  contents  of  the  vectors  in  (51),  consult  reference  [5].  To 
achieve  a  versatile  numerical  form  for  (47),  we  integrate  this  equation  over  a  finite  volume  cell  of 
volume  V. 


fffy§  dV  +  ffJyV-kv  =  0  (52) 

We  note  that  the  first  integral  is  time  independent,  so  for  this  term,  the  order  of  differentiation  and 
integration  may  be  interchanged.  To  the  second  integral,  the  divergence  theorem  is  applied;  hence, 


a 

dt 


Iffy  Q  +  ft(i/) 


F-dS  =  0 


(53) 


An  interesting  fact  concerning  (53)  is  that  the  second  integral  has  been  reduced  to  an  integral  over 
the  closed  surface  (denoted  d(V))  surrounding  the  finite  volume  cell.  The  leftmost  integral  in  (53) 
contains  the  volume  integral  of  the  vector  of  conserved  variables.  This  integral  can  be  envisioned 
as  the  volume  average  of  Q  multiplied  by  the  cell  volume.  With  this  assertion  in  mind,  let  us 
evaluate  (53)  for  cell  (i,j,  k).  We  obtain 


or 


V, 


i,j,k 


dQi,j,k 

dt 


"ft 


F  ■  dS 


dQi,j,k 

dt 


F  ■  dS 


(54) 


(55) 


This  procedure  has  reduced  the  system  of  governing  equations  to  a  system  of  ordinary  differential 
equations,  one  for  each  cell.  Recalling  that  the  finite  volume  cells  are  six-sided  (hexahedral),  the 
surface  integral  is  performed  over  six  planar  sides  indexed  by  s.  By  implementing  this  idea,  we 
have  that 


dQi,j,k 

dt 


•  nsAs 


(56) 


where  ns  is  the  outward  pointing  normal  vector  for  cell  side  ,v,  and  As  is  the  area  of  the  side.  The 

construct  Ffj  k is  the  flux  vector  evaluated  at  side  s  of  cell  (i ,  j,  k).  Equation  (56)  is  in  what  is 
referred  to  a  semi-discrete  form.  The  right  side  of  this  equation  is  discretized  pending  the 
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mechanics  of  the  space  scheme  while  the  left  side  is  an  ordinary  time  derivative.  Of  course,  the 
Vi j  k  are  the  cell  volumes.  The  left  side  still  remains  to  be  discretized. 

To  accurately  resolve  wave  motion  in  a  given  flow  field,  LESLIE3D  is  developed  with  an 
“explicit”  time  marching  algorithm.  That  is  to  say,  the  solution  at  time  step  N+l  relies  only  upon 
the  solution  computed  at  time  N.  This  property  of  the  solver  renders  highly  accurate  numerical 
solutions  for  waves  with  minimal  numerical  dissipation,  and  it  streamlines  the  structure  of  the 
code.  A  disadvantage  is  that  often  a  very  small  time  step  size  is  required  to  maintain  stability  for 
the  numerical  solution.  Of  course,  LESLIE3D  adjusts  its  time  step  automatically.  The  time 
marching  scheme  involved  is  a  relatively  simple  Runge-Kutta  method  devised  for  integrating 
ordinary  differential  equations.  [12]  If  we  regard  the  right  side  of  (54)  as  A  Q/ At,  then  this  scheme 
may  be  written  in  two  steps  as 


Qa  _  QU  +  AQn 

(57) 

Q n+1  =  +  (T  +  A(r] 

In  this  case,  Qa  is  the  intermediate  solution  between  times  n  and  n+l.  Given  the  small  time  steps 
required  for  stability,  this  time  stepping  scheme  is  very  accurate. 

The  larger  portion  of  the  work  associated  with  the  computation  of  a  LESLIE3D  solution 
lies  in  accomplishing  the  integration  in  space,  i.e.,  calculating  the  sum  in  equation  (56).  Years  ago, 
LESLIE3D  was  really  developed  with  the  use  of  second  and  fourth  order  MacCormack 
schemes. [5]  For  shock  free  flow  fields,  these  schemes  work  quite  well.  Possessing  very  little 
dissipation,  these  space  schemes  capture  turbulent  motions  in  the  flow  field  with  great  accuracy. 
Unfortunately,  when  LESLIE3D  was  adapted  for  shocked  flowfields,  the  MacCormack  scheme 
performed  poorly  (as  expected).  MacCormack  methods  are  basically  centered  finite  difference 
schemes  meaning  that  information  on  the  grid  is  weighted  symmetrically  in  all  directions.  Shocked 
flow  fields  contain  shock  waves,  and  in  the  vicinity  of  a  shock  wave  information  does  not  arrive 
from  all  directions.  Rather,  the  flow  of  information  has  a  preferential  direction.  In  fact,  it  comes 
from  upstream  of  the  shock  wave,  and  the  numerical  algorithms  to  capture  shock  waves  must  be 
written  accordingly.  The  associated  algorithms  are  called  upwind  algorithms  for  this  reason.  There 
are  many  upwind  algorithms  divided  into  two  major  subdivisions:  flux  vector  splitting  schemes 
and  flux  difference  splitting  schemes.  Generally,  for  problems  that  involve  more  complicated 
equations  of  state,  flux  difference  splitting  methods  are  preferred.  These  schemes  are  closely 
related  to  Godunov's  method,  but  where  Godonov's  scheme  entails  an  exact  solution  for  the 
Riemann  problem,  flux  difference  splitting  schemes  tend  to  use  approximate  Riemann 
solvers. [13, 14]  The  Riemann  problem  basically  addresses  the  flow  solution  that  occurs  at  a 
discontinuity.  A  discontinuity  can  be  thought  of  as  two  volumes  of  gas  separated  by  a  fictitious, 
massless  barrier  that  vanishes  at  the  start  of  the  problem.  On  one  side  of  the  barrier,  we  may  have 
high  pressure,  temperature  and  density  while  on  the  other  side,  these  properties  have  lower 
magnitudes.  We  apply  this  idea  at  every  cell  side  in  the  flow  field.  Approximate  Riemann  solvers, 
like  Roe's  solver  [15],  provide  an  exact  solution  for  this  problem  at  each  interface.  Roe's  method 
is  of  substantial  theoretical  interest,  but  it  difficult  to  apply  for  complicated  equations  of  state.  [16] 
Other  methods  that  do  not  so  heavily  rely  on  a  mechanical  characteristic  wave  decomposition  are 
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preferred  for  general  application. 

This  effort  has  encountered  great  success  with  variations  of  the  Harten,  Lax  and  van  Leer 
(HLL)  family  of  schemes.  Departing  from  the  strict  flux  difference  splitting  form  for  the  flux, 
i.e., 


dF  —  AdQ 


(58) 


where  A  is  a  flux  Jacobian  matrix,  the  Riemann  problem  is  sketched  as  an  entity  in  space-time 
coordinate  system.  Space,  in  this  regard,  is  cast  in  one  dimension.  For  a  simple  description  of  the 
HLL  solver,  suppose  a  discontinuity  is  initiated  at  x  =  0;  the  solution  is  formulated  by  one  left 
traveling  wave  and  one  right  traveling  wave.  [17]  The  x  direction,  as  it  were,  is  aligned 
perpendicular  to  the  discontinuity  (or  interface).  The  HLL  method  is  built  upon  an  integral  form 
of  the  Euler  equations.  Note  that  the  Euler  equations  are  considered  in  this  analysis  because  upwind 
methods  are  designed  for  advective  terms.  The  Navier-Stokes  equations  possess  the  same 
advective  terms.  The  integral  form  is  obtained  as  follows.  The  vector  form  (51)  of  our  system  of 
conservation  laws  can  be  written  aligned  perpendicular  to  an  interface  as 


dQ  dF 

dt  dx 


0 


(59) 


By  performing  a  cyclical  integral  around  the  rectangular  x-t  region,  we  have  that 

f  §  ~  I)  **  =  0  («» 

Hence, 

§(Qdx  -  Fdt )  =  0  (61) 

— >  — > 

The  multi-component,  non- geometric,  vectors  Qand  F  are  easily  evaluated  on  the  fixed  x  and  fixed 
t  rectangular  boundaries  of  this  region  resulting  in  a  system  of  equations  that  relate  their 
components  allowing  the  shock  properties  to  be  calculated.  To  the  left  of  the  left  traveling  wave, 

the  properties  are  uniform  with  values  QLwith  a  uniform  flux  vector  FL.  The  wave  is  assumed  to 

— )  — ^ 

propagate  with  speed  SL.  Corresponding  values  for  the  right  traveling  wave  are  QR,  FR and  SR, 
respectively.  The  region  in  between  the  left  and  right  traveling  waves  is  assumed  to  have  the 
uniform  properties  Q  and  flux  F* .  Equation  (61)  can  then  be  evaluated  on  the  perimeter  of  the 
rectangular  x-t  region  to  show  that 


_  fl-slql-(fr-srqr ) 
V  Sr-sl 


(62) 


as  the  properties  vector  in  the  region  between  the  waves.  After  applying  (61)  across  the  individual 
left  and  right  waves,  it  can  be  shown  that 
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(63) 


p*  _  srfl-slfr-slsr(qr-Ql) 

SR~SL 


is  the  flux  vector  for  the  region  between  the  waves.  The  upwind  flux  is  then  defined  at  the  interface 
(located  here  at  x  =  0,  but  at  the  cell  side  in  practice)  as 


(64) 


The  left  and  right  wave  speeds  must  still  be  estimated  before  this  method  is  utilized.  Commonly 
used  wave  speed  estimates  are  provided  by  the  work  of  Einfeldt.[18]  These  estimates  are 


SL  =  min [14  -  cL,  (V)  -  <c>] 


(65) 


SR  =  max [14  +  cR,  ( V )  +  <c>] 


(66) 


The  angle  brackets  used  in  (65)  and  (66)  are  indicative  of  the  use  of  a  Roe  average.  The  HLL 
method  is  not  difficult  to  implement,  but  it  does  have  one  disadvantage.  It  does  not  recognize  the 
presence  of  a  contact  discontinuity  between  the  waves.  Instead  the  solution  is  effectively  smeared 
across  this  region.  The  contact  discontinuity  is  added  back  into  the  Riemann  problem  solution  in 
the  HLLC  solver.  [19]  Due  to  the  additional  complexity  of  this  solver  and  in  keeping  with  the 
brevity  of  this  report,  HLLC  is  not  derived  here. 

There  are  many  important  aspects  of  upwind  solvers  that  cannot  be  included  in  this  brief 
report.  Among  these  are  interface  variable  reconstruction  schemes,  total  variation  diminishing 
(TVD)  schemes  (limiters)  and  others  that  allow  first  order  upwind  scheme  solutions  like  HLL  to 
be  extrapolated  to  higher  order.  Lor  those  readers  who  are  interested,  refer  to  References  [5],  [11] 
and  [13];  many  other  such  references  also  exist  in  the  public  domain. 

Another  core  capability  of  LESLIE3D's  algorithms  that  requires  mention  is  the  hybrid 
shock-turbulence  capture  scheme.  Earlier  in  the  report,  it  is  mentioned  that  the  algorithms  for 
capturing  shocks  clash  with  the  algorithms  for  capturing  turbulence  in  the  flow  field.  The  former 
adds  a  significant  amount  of  numerical  viscosity  while  the  latter  abhors  it.  LESLIE3D  possesses  a 
specially  modified  HLLC  shock-capturing  scheme  that  includes  some  of  the  features  of  the 
Piecewise  Parabolic  Method  (PPM)  [20]  to  create  an  approximate  Riemann  solver  that  has  very 
little  excess  numerical  dissipation.  Secondly,  LESLIE3D  adapts  a  special  switching  method  to 
only  apply  the  pure  upwind  solver  in  regions  where  discontinuities  are  detected.  In  smooth  regions 
of  flow,  the  MacCormack  scheme  is  employed  to  accurately  capture  turbulence.  The  upwind 
method  is  adjoined  carefully  to  the  MacCormack  scheme  through  a  special  shock  detection 
algorithm  to  ensure  that  the  appropriate  solver  algorithm  is  turned  on  at  the  appropriate  time. 
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Figure  1.  Spherical  grid  with  one  block  set  aside  for  illustrative  purposes 


Figure  2.  Perturbation  induced  along  an  arc  on  the  inner  surface  of  block  10 


3.5  Problem  Set-Up 

As  mentioned  near  the  beginning  of  the  this  report,  the  problem  at  hand  is  to  validate  post¬ 
processing  computer  codes  developed  in-house  against  physical  parameters  for  a  known  flow 
configuration.  Due  to  the  substantial  amount  of  available  drag  data,  flow  around  a  sphere  is 
selected  for  study,  and  the  key  parameter  for  comparison  is  the  average  drag  coefficient.  In  truth, 
we  expect  that  the  drag  coefficient  is  somewhat  time  dependent,  perhaps  with  the  some  periodic 
character  depending  on  the  Reynolds  number.  Naturally,  a  grid  must  be  generated  for  this 
configuration.  The  general  form  of  a  spherical  shell  is  chosen.  The  flow  field  grid  exists  between 
the  inner  and  outer  spherical  surfaces.  Figure  1  illustrates  this  concept.  The  grid  is  produced  by 
rotating  key  vectors  one  onto  another  with  the  origin  as  a  vertex.  In  particular,  an  octant  is  divided 
into  three  grid  volumes  by  selecting  a  central  vector  for  the  octant.  For  example,  the  first  octant  is 
delineated  by  the  vectors: 

^=(1,0,0);  y2=(0,l,0);  V3=(0,0,l)  (67) 

This  volume  can  be  divided  into  three  smaller  volumes  of  equal  size  and  shape  by  adding  the 
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vector 


/ 


V  = 


\ 


1/  1/  1 

\7S'  7S’7Sj 


(68) 


as  the  central  vector.  Accomplishing  this  procedure  in  each  of  the  eight  octants  between  the  inner 
and  outer  spherical  shells  leads  to  the  formation  of  24  volumes  or  blocks.  To  facilitate  viscous 
spacing,  the  grid  is  stretched  in  the  radial  direction  equivalently  in  each  block.  A  simple  finite 
geometric  series  is  used  as  the  stretching  function.  The  spherical  shell  grid  has  inner  and  outer 
radii  of  0.1  meter  and  1  meter,  respectively.  For  this  test  case,  the  minimum  spacing  at  the  sphere's 
surface  is  set  at  10"4  meter.  To  provide  for  LESLIE3D's  multi-block  data  communication,  two 
ghost  surfaces  are  installed  adjacent  to  each  of  the  six  boundary  surfaces  for  each  block.  For  those 
surfaces  that  communicate  to  an  adjacent  block,  the  ghost  surfaces  must  precisely  overlay  on  the 
adjacent  block's  grid  surfaces.  If  there  is  a  mismatch  in  the  ghost  point  locations,  errors  quickly 
arise  in  LESLIE3D's  derivative  computations.  A  final  issue  that  sometimes  arises  for  this  type  of 
problem  is  that  of  flow  symmetry.  The  geometry  of  this  grid  configuration  is  highly  symmetric.  In 
some  cases,  this  symmetry  can  be  inherited  by  the  numerical  solution  creating  an  unlikely  physical 
realization  of  the  flow  field  around  the  sphere.  Two  measures  are  implemented  to  prevent  this 
problem  from  occurring.  First,  the  grid  is  altered  at  the  sphere's  surface  in  block  10  on  the  -x  side 
of  the  sphere.  A  perturbation  is  made  along  an  arc  that  extends  across  the  block  face;  at  the  block 
boundary,  the  perturbation  is  zero  while  at  the  center  of  the  arc,  the  perturbation  is  a  maximum. 
The  shape  of  the  perturbation  is  that  of  an  arc-shaped  indentation  on  the  sphere.  A  picture  of  this 
perturbation  is  shown  in  Figure  2. 

Another  important  aspect  of  computing  a  flow  solution  for  the  sphere  is  that  of  boundary 
conditions.  As  one  may  expect,  boundary  conditions  have  to  be  established  for  this  problem  in 
each  block.  The  lateral  surfaces,  those  with  constant  index  coordinates,  (i  =  1,  IMAX  or  k  =  1, 
KMAX)  are  communicating  interfaces  that  exchange  data  with  neighboring  blocks.  This  type  of 
condition  is  not  a  true  boundary  condition.  The  inner  surface  boundary  is  quite  simple.  A  viscous 
(no-slip)  boundary  condition  is  enforced  for  the  inner  (j  =  1)  surface  for  each  block.  The  outer 
boundary  is  a  bit  more  complicated.  The  yz-planc  separates  the  regions  of  inflow  from  the  regions 
of  outflow.  Given  the  order  in  which  block  grids  are  computed,  the  outer  (j  =  JMAX)  boundary 
for  blocks  1  through  12  is  designated  as  an  inflow.  The  JMAX  boundary  for  the  remaining  blocks 
13  through  24  is  designated  as  a  supersonic  outflow.  FESFIE3D  allows  the  user  to  specify  the  type 
of  inflow  and  outflow  as  either  characteristic  variable  for  subsonic  flow  or  interpolated  for 
supersonic  flow.  Although  the  boundary  conditions  are  specified  in  an  input  file  by  the  user,  these 
conditions  are  read  by  the  pre-processor  (pre_les3d)  and  formally  set  in  LESLIE3D's  processor 
communication  files. 


LESLIE3D  also  requires  that  the  user  specify  initial  conditions.  For  the  purposes  of  this 
work,  initial  conditions  are  specified  in  an  initial  restart  file  REST_000.DATA  created  in  the  data 
directory  for  each  block.  Pressure,  temperature,  velocity  components,  subgrid  kinetic  energy  and 
initial  chemical  species  mass  fractions  must  be  set  by  the  user.  For  the  problems  solved  here, 
uniform  flow  conditions  are  set  for  the  entire  flow  field  without  regard  for  the  presence  of  the 
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Figure  3.  Grid  block  arrangement  for  the  simulation  of  subsonic  flow  outside  of  the  sphere 


sphere's  surface.  Starting  the  numerical  solution  in  this  way  creates  a  "slug  flow"  condition  that  is 
non-physical.  Solver  time  is  required  for  LESLIE3D  to  advect  this  non-physical  slug  of  flow  from 
the  solution.  For  refined  grids,  this  process  can  be  time  consuming.  An  approach  used  in  this  case 
is  to  first  run  the  problem  on  a  uniform  (equally  spaced)  radial  grid.  This  process  allows  larger 
time  steps  and  can  more  quickly  dispense  with  the  slug  flow.  At  this  point,  the  numerical  solution 
is  remapped  onto  a  refined  grid,  and  LESLIE3D  is  restarted.  The  solution  is,  of  course,  incorrect 
near  the  body  surface  due  to  the  lack  of  refinement,  so  LESLIE3D  is  run  forward  with  the  LES 
algorithms  activated  until  the  flow  near  the  surface  becomes  stationary.  Stationarity  can  be  verified 
by  examining  time  histories  for  flow  variables  say,  vorticity  magnitude  and  subgrid  kinetic  energy. 
A  series  of  restart  files  taken  after  the  flow  solution  has  become  stationary  may  be  used  to  calculate 
the  force  time  histories,  e.g.  the  drag  force  in  time. 

The  spherical  grid  presented  above  functions  well  for  supersonic  flow  problems  but  not  for 
subsonic  flow.  At  subsonic  speeds,  the  treatment  of  inflow  and  outflow  boundaries  requires  special 
attention.  At  these  boundaries,  acoustic  waves  can  enter  or  exit  the  computational  domain.  As  a 
result,  these  waves  can  alter  the  numerical  solution.  A  better  grid  is  generated  for  this  problem  by 
projecting  the  outer  surfaces  of  the  component  spherical  blocks  onto  flat  planes.  Rectangular 
inflow  and  outflow  blocks,  numbered  25  through  32,  are  added  to  the  grid  to  implement 
characteristic  variable  boundary  conditions.  This  grid  arrangement  is  shown  in  Figure  3.  The  blue 
block  shown  in  Figure  3  connects  the  sphere  to  one  of  brick  shaped  inflow  blocks.  Similar  blocks 
connect  to  the  brick-shaped  outflow  blocks.  The  red  block  extends  laterally  to  the  outer  boundary 
which  is  coded  as  a  slip  surface.  LESLIE3D’s  characteristic  variable  inflow/outflow  boundary 
conditions  require  grid  blocks  with  a  higher  degree  of  uniformity  and  more  importantly,  coincident 
indicial  axis  directions.  The  revised  grid  design  satisfies  these  requirements. 
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4.0  RESULTS 


Figure  4.  Slice  pressure  plot  on  the  xy-plane  of  the  sphere  flow  field  at  Mach  2,  units  in  Pascals 
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Figure  5.  Contoured  slice  pressure  plot  on  the  *y-plane  of  the  sphere  flow  field  at  Mach  2,  units  in  Pascals 
4.1  Supersonic  Test  Case  -  Flow  Field  Results 

The  first  test  case  considered  assumes  that  the  sphere  is  traveling  at  Mach  2,  altitude  1000 
feet.  The  altitude  data  allows  the  freestream  pressure  and  density  to  be  set  based  upon  atmospheric 
tables.  The  Reynolds  number  is  computed  as  about  8.1  million  based  upon  the  sphere's  diameter 
of  0.2  meters.  At  this  altitude,  the  freestream  velocity  is  667.7  meters  per  second.  Figure  4  contains 
a  slice  plot  of  the  pressure  field  taken  along  the  xy-plane.  The  curved  shock  wave  is  evident  as  is 
the  wake  region  even  though  the  grid  is  not  zoned  well  in  the  far  field.  Rather  the  best  zoning  is 
adjacent  to  the  body  to  support  better  drag  calculations.  Figure  5  contains  a  contoured  slice  plot  of 
pressure  taken  from  the  flow  field  along  the  xy-plane  at  12.5  ms  solution  time.  In  this  case,  the 
region  near  the  sphere's  surface  is  magnified  for  better  viewing.  The  stagnation  region  is  clearly 
visible,  and  the  wake  region  is  characterized  by  asymmetric  contouring,  just  as  expected. 
Temperature  plots  show  similar  behavior.  A  slice  plot  along  the  xy-plane  is  presented  in  Figure  6 
with  absolute  temperature  indicated  in  Kelvin.  The  high  temperature  stagnation  region  is  evident, 
but  more  interestingly,  the  gas  heated  in  the  boundary  layer  washes  off  from  the  sphere's  surface 
and  advects  into  the  wake.  This  phenomenon  is  also  visible  in  Figure  7,  a  close  in  view  of  the  same 
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slice. 


Figure  7.  Contoured  slice  temperature  plot  of  the  sphere  flow  field  at  Mach  2,  units  in  Kelvin 


Figure  8.  Slice  plot  of  u  velocity  with  contours  of  velocity  magnitude  at  13.61  ms;  units  are  m/s 
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Figure  9.  Slice  plot  of  v  velocity  with  contours  of  velocity  magnitude  at  13.61  ms;  units  are  m/s 


Figure  10.  Slice  plot  of  w  velocity  with  contours  of  velocity  magnitude  at  13.61  ms;  units  are  m/s 
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Figure  11.  Time  history  plots  for  x-component  of  vorticity,  units  s'1 
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Figure  12.  Time  history  plots  for  y-component  of  vorticity,  units  s'1 


2.5x10  6 
2x10  6 
1.5x10  6 
IxlO6 
500000 
0 

-500000 
-IxlO6 
-1.5x106 
-2x1 06 
-2.5x106 

0.0115  0.012  0.0125  0.013  0.0135  0.014  0.0145  0.015  0.0155  0.016 

Time  (s) 

B01  -  B05  B09  B13  B17  -  B21 

B02  -  B06  B10  B14  -  B18  -  B22  - 

B03  B0  7  -  B11  B15  -  B19  B23  - 

B04  -  B08  B12  B16  -  B20  -  B24  - 

Figure  13.  Time  history  plots  for  z-component  of  vorticity,  units  s'1 

Figure  8  is  an  interesting  slice  plot  of  u  velocity,  the  component  of  flow  in  the  stream  (or 
x )  direction.  Contour  lines  of  velocity  magnitude  are  applied  to  the  plot.  Flow  separation  is  evident 
since  gas  is  observed  propagating  upstream  in  the  region  near  the  back  of  the  sphere.  A 
recirculation  region  is  likely  to  exist  in  the  wake.  This  idea  is  reinforced  by  Figure  9,  a  similarly 
designed  plot  of  v  velocity,  the  vertical  component.  In  this  case,  we  can  see  masses  of  fluid 
alternatively  rising  and  falling  in  the  wake  indicative  of  a  recirculating  flow,  at  least  intermittently. 
This  idea  is  also  borne  out  by  Figure  10,  a  plot  of  w  velocity,  the  component  in  the  z  direction.  In 
this  case,  bright  red  color  indicates  flow  out  of  the  page  while  blue  is  in  the  opposite  direction. 
From  this  evidence,  it  is  obvious  that  this  flow  field  is  fully  three-dimensional. 
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Figure  14.  Time  history  plots  for  vorticity  magnitude,  units  s'1 
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Figure  15.  Time  history  plots  for  subgrid  kinetic  energy,  units  m2/s2 

Prior  to  post-processing  solver  output  data  for  drag  information,  it  is  necessary  to  show 
that  the  flow  field  is  stationary.  That  is  to  say,  its  statistical  quantities  are  unchanging.  We  can 
illustrate  this  idea  by  the  use  of  time  histories  of  fluctuating  quantities  like  vorticity  or  subgrid 
kinetic  energy.  Also,  we  examine  how  these  properties  behave  on  the  sphere's  surface  in  time.  If 
these  properties  exhibit  periodic  behavior  in  time,  then  the  flow  is  stationary.  To  examine  the 
approach  to  stationarity,  we  consider  flx,  fly ,  flz  and  fl,  the  x,  y  and  z  components  of  vorticity 
along  with  the  vorticity  magnitude.  Also,  the  subgrid  kinetic  energy  ksas  is  studied.  These 
properties  are  selected  because  of  their  strong  roles  in  the  evolution  of  turbulence.  Time  histories 
are  sampled  on  the  sphere's  surface  with  exactly  one  sampling  point  per  block.  In  most  cases,  the 
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Figure  16.  Drag  coefficient  time  history  for  the  sphere  at  Mach  2.0 


sampling  point  is  located  at  the  center  of  the  block's  inner  surface.  The  two  exceptions  are  blocks 
12  and  18,  where  the  sample  point  is  located  near  the  block  boundary.  Time  histories  for  fix,  fly, 
flz  and  (2  are  shown  in  Figures  11,  12,  13  and  14,  respectively.  Time  evolution  of  ksgs  is  shown  in 
Figure  15.  In  each  of  these  figures,  we  observe  that  envelopes  form  around  the  time  history  traces 
for  almost  all  of  the  sample  points.  The  exceptions  are  the  x-component  of  vorticity  and  with 
subgrid  scale  kinetic  energy.  Traces  for  fix  depart  from  the  envelope  for  sample  points  23  and  24. 
These  points  are  located  on  the  trailing  section  of  the  sphere  in  the  +y,+z  octant.  The  excursion  in 
occurs  at  sample  point  14  on  the  trailing  section  of  the  sphere  in  the  +x,-y,-z  octant.  Both  of 
these  regions  are  subject  to  vortex  shedding  and  recirculation,  so  it  is  expected  that  both  of  these 
properties  rise  and  fall  in  time.  Histories  formed  at  the  remainder  of  the  sample  points  exhibit 
nearly  periodic  behavior.  Based  upon  these  plots,  the  flow  field  is  sufficiently  stationary  to  admit 
force  calculations.  Force  components  (particularly  drag,  the  force  exerted  in  the  +x  direction)  are 
obtained  by  post-processing  a  series  of  solution  files.  For  the  supersonic  solution,  a  time  history 
for  the  drag  coefficient  is  shown  in  Figure  16.  The  archived  drag  coefficient  for  this  configuration 
is  essentially  1.0,  so  the  numerical  solution  agrees  well  with  this  datum.[22] 

4.2  Subsonic  Test  Case  -  Flow  Field  Results 

The  previous  test  case  examines  the  flow  field  around  the  sphere  at  high  Reynolds  number 
where  the  drag  curve  levels  out.  In  itself,  that  test  case  is  very  important  because  high  Reynolds 
number  situations  are  widely  applicable.  Still,  it  is  important  to  verify  the  performance  of  our 
computer  codes  at  low  Reynolds  number  where  the  drag  curve  exhibits  changes.  At  a  Reynolds 
number  of  approximately  400,000,  the  drag  curve  dips  sharply  due  to  boundary  layer  separation.  [3] 
It  is  desirable  to  verify  our  computer  codes  for  this  case.  The  flow  geometry  and  grids  remain  the 
same  as  in  the  preceding  case.  The  freestream  flow  conditions  do  change.  This  problem  is 
conducted  at  an  altitude  of  1,000  feet  and  Mach  0.1,  a  low  speed  flow  condition  corresponding  to 
about  33.3  m/s.  A  major  change  in  problem  set-up  involves  LESLIE3D's  boundary  conditions.  At 
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Figure  17.  Slice  plot  of  the  pressure  field  for  the  subsonic  solution  at  43.68  ms 


Figure  18.  Slice  plot  of  the  temperature  field  for  the  subsonic  solution  at  43.68  ms 


this  flow  speed,  acoustic  waves  can  travel  in  both  directions,  with  and  against  the  flow.  Physically, 
this  means  that  an  acoustic  signal  can  travel  from  the  sphere  to  the  numerical  inflow  boundary 
altering  the  flow  at  that  location.  This  situation  cannot  happen  in  supersonic  flow,  so  different 
algorithms  must  be  used  to  account  for  this  type  of  wave  propagation.  To  properly  capture  physics 
at  the  inflow  boundary,  characteristic  variable  boundary  conditions  are  employed.  [21]  This  type 
of  algorithm  uses  entropy,  vorticity  and  acoustic  wave  information  from  the  different  flow 
directions  to  solve  simultaneously  for  the  inflow  conditions.  The  outflow  boundary  is  also 
structured  in  a  similar  manner.  Characteristic  boundary  conditions  are  employed  at  the  outflow  to 
set  properties  at  the  outflow  plane.  Sample  plots  of  the  flow  field  are  shown  below  at  solution  time 
43.68  ms.  Figure  17  contains  a  slice  plot  of  pressure  while  Figure  1 8  contains  a  plot  of  temperature. 
MacCormack’s  space  scheme  is  utilized  for  computing  the  subsonic  solution  since  the  flow  field 
is  turbulent  but  completely  free  of  shock  waves.  Vortex  shedding  is  evident  in  both  Figures  17  and 
18;  the  downstream  direction  is  to  the  right  in  both  figures.  Note  the  high  pressure  stagnation  point 
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Figure  19.  Drag  coefficient  time  history  for  the  sphere  at  Mach  0.1 


existing  on  the  windward  side  of  the  sphere  in  Figure  17.  The  elevated  temperature  generated  in 
the  boundary  layer  is  shed  into  the  flow  field  at  the  back  of  the  sphere.  A  high  temperature  locus 
corresponding  to  this  region  is  visible  in  Figure  18.  Vortices  shed  somewhat  irregularly  from  the 
sphere  creating  oscillatory  forces,  and  the  drag  force  may  be  extracted  by  post-processing. 
Accordingly,  a  time  history  for  the  drag  coefficient  is  shown  in  Figure  19.  The  drag  coefficient  is 
unsteady,  but  its  average  is  calculated  as  0.464.  Note  that  the  early  time  values  (less  than  0.02  s) 
are  taken  before  the  solution  has  become  stationary.  The  average  Cd  for  this  Reynolds  number 
agrees  quite  well  with  the  archived  coefficient  plot  recorded  in  Schlicting.[3] 
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5.0  CONCLUSIONS 


This  report  focuses  on  a  validation  study  conducted  for  the  Large  Eddy  Simulation  with 
Linear  Eddy  modeling  in  3  Dimensions  (LESLIE3D)  multiphase  physics  computer  program.  This 
computer  program  has  application  for  air-to-air  missile  aerodynamics,  so  it  is  important  to 
demonstrate  user  capabilities  for  LESLIE3D  and  also  to  ensure  that  solution  data  can  be  accurately 
post-processed  to  extract  aerodynamics  data.  Aerodynamic  drag  is  of  principal  interest  in  this 
report.  The  validation  problem  selected  for  analysis  is  the  turbulent,  multi-species  flow  around  a 
sphere,  a  typical  bluff  body  often  used  for  validation.  The  flow  field  is  computed  first  at  Mach  2.0, 
Reynolds  number  8  million  and  then  at  Mach  0.1,  Reynolds  number  400,000.  The  first  problem 
involves  a  shocked  flow  field  while  the  latter  flow  field  is  smooth.  Both  flow  fields  are  turbulent 
and  set  at  an  altitude  of  1,000  feet.  In  each  case,  a  significant  amount  of  work  is  invested  in  the 
generation  of  grids  for  these  two  flow  fields.  In  each  case,  the  appropriate  numerical  boundary 
conditions  are  applied  for  the  two  different  flight  conditions.  Drag  coefficients  are  computed  from 
the  numerical  solution  for  these  problems.  For  the  supersonic  problem,  the  average  drag  coefficient 
is  computed  as  0.973,  and  for  the  subsonic  case,  the  average  drag  coefficient  is  calculated  as  0.464. 
These  results  compare  favorably  with  archived  data.  Accordingly,  these  results  serve  to  validate 
the  methodology  for  problem  set-up  employed  in  this  work.  Our  post-processing  algorithms  are 
validated  simultaneously.  Validation  problems  recommended  for  future  accomplishment  may 
involve  the  oblate  spheroids  or  hemisphere-cylinder  configurations.  These  problems  may  be 
solved  for  hypersonic  flight  conditions  with  atmospheric  chemical  reactions. 
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LIST  OF  ACRONYMS,  ABBREVIATIONS,  AND  SYMBOLS 


This  section  contains  brief  definitions  of  various  terms  and  acronyms  used  throughout  this 
document.  Only  terms  and  acronyms  whose  definitions  are  considered  uncommon  are  included. 


HLLC . Harten,  Lax  and  van  Leer  Contact  preserving 

LDKM . Locally  Dynamic  subgrid  Kinetic  Energy  Model 

LES . Large  Eddy  Simulation 

LESLIE3D . Large  Eddy  Simulation  with  Linear  Eddy  modeling  in  3  Dimensions 

MUSCL . Monotone  Upstream  centered  Schemes  for  Conservation  Laws 


Cd . Drag  coefficient 

Cp . Constant  specific  heat  at  constant  pressure 

Cv . Constant  volume  specific  heat  capacity 

Dk . Species  diffusion  coefficient 

E . Energy  per  unit  volume 

e . Internal  energy  per  unit  mass 

h . Sensible  enthalpy 

MWn . Molecular  weight  of  species  n 

P . Pressure 

qi . Heat  flux  vector  component 

R . Species  gas  constant 

Runiv . Universal  gas  constant 

Re . Reynolds  number 

sgs . Subgrid  scale 

T . Temperature 

m . Cartesian  velocity  component 

V  . Velocity  magnitude 

V  . Velocity  vector 

Vik . Diffusion  velocity 

Yk . Species  mass  fraction 


Xi.. 

X  . 

K  . 

5  . 
p- 

Tij 

p- 

(Ok 


Cartesian  space  coordinate 

Bulk  viscosity 

Heat  transfer  coefficient 

Delta  function 

Density 

Shear  stress  tensor 
Dynamic  viscosity 
Species  production  rate 
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